/*-------------------------------------------------------------------------
 *
 * geo-ops.c--
 *    2D geometric operations
 *
 * Copyright (c) 1994, Regents of the University of California
 *
 *
 * IDENTIFICATION
 *    $Header: /usr/local/cvsroot/postgres95/src/backend/utils/adt/geo-ops.c,v 1.3 1996/11/03 06:53:06 scrappy Exp $
 *
 *-------------------------------------------------------------------------
 */
#include <math.h>
#include <float.h>	/* faked on sunos */
#include <stdio.h>	/* for sprintf proto, etc. */
#include <string.h>

#include "postgres.h"

#include "utils/geo-decls.h"
#include "utils/palloc.h"

#define LDELIM		'('
#define RDELIM		')'
#define	DELIM		','
#define BOXNARGS	4
#define	LSEGNARGS	4
#define	POINTNARGS	2

/***********************************************************************
 **
 ** 	Routines for two-dimensional boxes.
 **
 ***********************************************************************/

/*----------------------------------------------------------
 * Formatting and conversion routines.
 *---------------------------------------------------------*/

/*	box_in	-	convert a string to internal form.
 *
 *	str:	input string "(f8, f8, f8, f8)"
 */
BOX *box_in(char *str)
{
    double	tmp;
    char	*p, *coord[BOXNARGS];
    int	i;
    BOX	*result;
    
    if (str == NULL)
	elog (WARN," Bad (null) box external representation");
    
    if ((p = (char *)strchr(str, LDELIM)) == (char *)NULL)
	elog (WARN, "Bad box external representation '%s'",str);
    for (i = 0, p = str; *p && i < BOXNARGS && *p != RDELIM; p++)
	if (*p == DELIM || (*p == LDELIM && !i))
	    coord[i++] = p + 1;
    if (i < BOXNARGS - 1)
	elog (WARN, "Bad box external representation '%s'", str);
    result = PALLOCTYPE(BOX);
    result->xh = atof(coord[0]);
    result->yh = atof(coord[1]);
    result->xl = atof(coord[2]);
    result->yl = atof(coord[3]);
    if (result->xh < result->xl) {
	tmp = result->xh;
	result->xh = result->xl;
	result->xl = tmp;
    }
    if (result->yh < result->yl) {
	tmp = result->yh;
	result->yh = result->yl;
	result->yl = tmp;
    }
    
    return(result);
}

/*	box_out	-	convert a box to external form.
 */
char *box_out(BOX *box)
{
    char	*result;
    
    if (box == NULL)
	return(NULL);
    result = (char *)PALLOC(80);
    (void) sprintf(result, "(%G,%G,%G,%G)",
		   box->xh, box->yh, box->xl, box->yl);
    
    return(result);
}


/*	box_construct	-	fill in a new box.
 */
BOX *box_construct(double x1, double x2, double y1, double y2)
{
    BOX	*result;
    
    result = PALLOCTYPE(BOX);
    return( box_fill(result, x1, x2, y1, y2) );
}


/*	box_fill	-	fill in a static box
 */
BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2)
{
    double	tmp;
    
    result->xh = x1;
    result->xl = x2;
    result->yh = y1;
    result->yl = y2;
    if (result->xh < result->xl) {
	tmp = result->xh;
	result->xh = result->xl;
	result->xl = tmp;
    }
    if (result->yh < result->yl) {
	tmp = result->yh;
	result->yh = result->yl;
	result->yl = tmp;
    }
    
    return(result);
}


/*	box_copy	-	copy a box
 */
BOX *box_copy(BOX *box)
{
    BOX	*result;
    
    result = PALLOCTYPE(BOX);
    memmove((char *) result, (char *) box, sizeof(BOX));
    
    return(result);
}


/*----------------------------------------------------------
 *  Relational operators for BOXes.
 *	<, >, <=, >=, and == are based on box area.
 *---------------------------------------------------------*/

/*	box_same	-	are two boxes identical?
 */
long box_same(BOX *box1, BOX *box2)
{
    return((box1->xh == box2->xh && box1->xl == box2->xl) &&
	   (box1->yh == box2->yh && box1->yl == box2->yl));
}

/*	box_overlap	-	does box1 overlap box2?
 */
long box_overlap(BOX *box1, BOX *box2)
{
    return(((box1->xh >= box2->xh && box1->xl <= box2->xh) ||
	    (box2->xh >= box1->xh && box2->xl <= box1->xh)) &&
	   ((box1->yh >= box2->yh && box1->yl <= box2->yh) ||
	    (box2->yh >= box1->yh && box2->yl <= box1->yh)) );
}

/*	box_overleft	-	is the right edge of box1 to the left of
 *				the right edge of box2?
 *
 *	This is "less than or equal" for the end of a time range,
 *	when time ranges are stored as rectangles.
 */
long box_overleft(BOX *box1, BOX *box2)
{
    return(box1->xh <= box2->xh);
}

/*	box_left	-	is box1 strictly left of box2?
 */
long box_left(BOX *box1, BOX *box2)
{
    return(box1->xh < box2->xl);
}

/*	box_right	-	is box1 strictly right of box2?
 */
long box_right(BOX *box1, BOX *box2)
{
    return(box1->xl > box2->xh);
}

/*	box_overright	-	is the left edge of box1 to the right of
 *				the left edge of box2?
 *
 *	This is "greater than or equal" for time ranges, when time ranges
 *	are stored as rectangles.
 */
long box_overright(BOX *box1, BOX *box2)
{
    return(box1->xl >= box2->xl);
}

/*	box_contained	-	is box1 contained by box2?
 */
long box_contained(BOX *box1, BOX *box2)
{
    return((box1->xh <= box2->xh && box1->xl >= box2->xl &&
	    box1->yh <= box2->yh && box1->yl >= box2->yl));
}

/*	box_contain	-	does box1 contain box2?
 */
long box_contain(BOX *box1, BOX *box2)
{
    return((box1->xh >= box2->xh && box1->xl <= box2->xl &&
	    box1->yh >= box2->yh && box1->yl <= box2->yl));
}


/*	box_positionop	-
 *		is box1 entirely {above, below } box2?
 */
long box_below(BOX *box1, BOX *box2)
{
    return( box1->yh <= box2->yl );
}

long box_above(BOX *box1, BOX *box2)
{
    return( box1->yl >= box2->yh );
}


/*	box_relop	-	is area(box1) relop area(box2), within
 *			  	our accuracy constraint?
 */
long box_lt(BOX *box1, BOX *box2)
{
    return( FPlt(box_ar(box1), box_ar(box2)) );
}

long box_gt(BOX *box1, BOX *box2)
{
    return( FPgt(box_ar(box1), box_ar(box2)) );
}

long box_eq(BOX *box1, BOX *box2)
{
    return( FPeq(box_ar(box1), box_ar(box2)) );
}

long box_le(BOX	*box1, BOX *box2)
{
    return( FPle(box_ar(box1), box_ar(box2)) );
}

long box_ge(BOX	*box1, BOX *box2)
{
    return( FPge(box_ar(box1), box_ar(box2)) );
}


/*----------------------------------------------------------
 *  "Arithmetic" operators on boxes.
 *	box_foo	returns foo as an object (pointer) that
 can be passed between languages.
 *	box_xx	is an internal routine which returns the
 *		actual value (and cannot be handed back to
 *		LISP).
 *---------------------------------------------------------*/

/*	box_area	-	returns the area of the box.
 */
double *box_area(BOX *box)
{
    double *result;
    
    result = PALLOCTYPE(double);
    *result = box_ln(box) * box_ht(box);
    
    return(result);
}


/*	box_length	-	returns the length of the box 
 *				  (horizontal magnitude).
 */
double *box_length(BOX *box)
{
    double	*result;
    
    result = PALLOCTYPE(double);
    *result = box->xh - box->xl;
    
    return(result);
}


/*	box_height	-	returns the height of the box 
 *				  (vertical magnitude).
 */
double *box_height(BOX *box)
{
    double	*result;
    
    result = PALLOCTYPE(double);
    *result = box->yh - box->yl;
    
    return(result);
}


/*	box_distance	-	returns the distance between the
 *				  center points of two boxes.
 */
double *box_distance(BOX *box1, BOX *box2)
{
    double	*result;
    Point	*box_center(), *a, *b;
    
    result = PALLOCTYPE(double);
    a = box_center(box1);
    b = box_center(box2);
    *result = HYPOT(a->x - b->x, a->y - b->y);
    
    PFREE(a);
    PFREE(b);
    return(result);
}


/*	box_center	-	returns the center point of the box.
 */
Point *box_center(BOX *box)
{
    Point	*result;
    
    result = PALLOCTYPE(Point);
    result->x = (box->xh + box->xl) / 2.0;
    result->y = (box->yh + box->yl) / 2.0;
    
    return(result);
}


/*	box_ar	-	returns the area of the box.
 */
double box_ar(BOX *box)
{
    return( box_ln(box) * box_ht(box) );
}


/*	box_ln	-	returns the length of the box 
 *				  (horizontal magnitude).
 */
double box_ln(BOX *box)
{
    return( box->xh - box->xl );
}


/*	box_ht	-	returns the height of the box 
 *				  (vertical magnitude).
 */
double box_ht(BOX *box)
{
    return( box->yh - box->yl );
}


/*	box_dt	-	returns the distance between the
 *			  center points of two boxes.
 */
double box_dt(BOX *box1, BOX *box2)
{
    double	result;
    Point	*box_center(),
    *a, *b;
    
    a = box_center(box1);
    b = box_center(box2);
    result = HYPOT(a->x - b->x, a->y - b->y);
    
    PFREE(a);
    PFREE(b);
    return(result);
}

/*----------------------------------------------------------
 *  Funky operations.
 *---------------------------------------------------------*/

/*	box_intersect	-
 *		returns the overlapping portion of two boxes,
 *		  or NULL if they do not intersect.
 */
BOX *box_intersect(BOX	*box1, BOX *box2)
{
    BOX	*result;
    long	box_overlap();
    
    if (! box_overlap(box1,box2))
	return(NULL);
    result = PALLOCTYPE(BOX);
    result->xh = Min(box1->xh, box2->xh);
    result->xl = Max(box1->xl, box2->xl);
    result->yh = Min(box1->yh, box2->yh);
    result->yl = Max(box1->yl, box2->yl);
    
    return(result);
}


/*	box_diagonal	-	
 *		returns a line segment which happens to be the
 *		  positive-slope diagonal of "box".
 *		provided, of course, we have LSEGs.
 */
LSEG *box_diagonal(BOX *box)
{
    Point	p1, p2;
    
    p1.x = box->xh;
    p1.y = box->yh;
    p2.x = box->xl;
    p2.y = box->yl;
    return( lseg_construct( &p1, &p2 ) );
    
}

/***********************************************************************
 **
 ** 	Routines for 2D lines.
 **		Lines are not intended to be used as ADTs per se,
 **		but their ops are useful tools for other ADT ops.  Thus,
 **		there are few relops.
 **
 ***********************************************************************/

/*----------------------------------------------------------
 *  Conversion routines from one line formula to internal.
 *	Internal form:	Ax+By+C=0
 *---------------------------------------------------------*/

LINE *				/* point-slope */
line_construct_pm(Point *pt, double m)
{
    LINE	*result;
    
    result = PALLOCTYPE(LINE);
    /* use "mx - y + yinter = 0" */
    result->A = m;
    result->B = -1.0;
    result->C = pt->y - m * pt->x;
    return(result);
}


LINE *				/* two points */
line_construct_pp(Point *pt1, Point *pt2)
{
    LINE	*result;
    
    result = PALLOCTYPE(LINE);
    if (FPeq(pt1->x, pt2->x)) {		/* vertical */
	/* use "x = C" */
	result->m = 0.0;
	result->A = -1.0;
	result->B = 0.0;
	result->C = pt1->x;
    } else {
	/* use "mx - y + yinter = 0" */
	result->m = (pt1->y - pt2->y) / (pt1->x - pt2->x);
	result->A = result->m;
	result->B = -1.0;
	result->C = pt1->y - result->m * pt1->x;
    }
    return(result);
}


/*----------------------------------------------------------
 *  Relative position routines.
 *---------------------------------------------------------*/

long line_intersect(LINE *l1, LINE *l2)
{
    return( ! line_parallel(l1, l2) );
}

long line_parallel(LINE *l1, LINE *l2)
{
    return( FPeq(l1->m, l2->m) );
}

long line_perp(LINE *l1, LINE *l2)
{
    if (l1->m)
	return( FPeq(l2->m / l1->m, -1.0) );
    else if (l2->m)
	return( FPeq(l1->m / l2->m, -1.0) );
    return(1);	/* both 0.0 */
}

long line_vertical(LINE *line)
{
    return( FPeq(line->A, -1.0) && FPzero(line->B) );
}

long line_horizontal(LINE *line)
{
    return( FPzero(line->m) );
}


long line_eq(LINE *l1, LINE *l2)
{
    double	k;
    
    if (! FPzero(l2->A))
	k = l1->A / l2->A;
    else if (! FPzero(l2->B))
	k = l1->B / l2->B;
    else if (! FPzero(l2->C))
	k = l1->C / l2->C;
    else
	k = 1.0;
    return( FPeq(l1->A, k * l2->A) &&
	   FPeq(l1->B, k * l2->B) &&
	   FPeq(l1->C, k * l2->C) );
}


/*----------------------------------------------------------
 *  Line arithmetic routines.
 *---------------------------------------------------------*/

double * 		/* distance between l1, l2 */
line_distance(LINE *l1, LINE *l2)
{
    double	*result;
    Point	*tmp;
    
    result = PALLOCTYPE(double);
    if (line_intersect(l1, l2)) {
	*result = 0.0;
	return(result);
    }
    if (line_vertical(l1))
	*result = fabs(l1->C - l2->C);
    else {
	tmp = point_construct(0.0, l1->C);
	result = dist_pl(tmp, l2);
	PFREE(tmp);
    }
    return(result);
}

Point *			/* point where l1, l2 intersect (if any) */
line_interpt(LINE *l1, LINE *l2)
{
    Point	*result;
    double	x;
    
    if (line_parallel(l1, l2))
	return(NULL);
    if (line_vertical(l1))
	result = point_construct(l2->m * l1->C + l2->C, l1->C);
    else if (line_vertical(l2))
	result = point_construct(l1->m * l2->C + l1->C, l2->C);
    else {
	x = (l1->C - l2->C) / (l2->A - l1->A);
	result = point_construct(x, l1->m * x + l1->C);
    }
    return(result);
}

/***********************************************************************
 **
 ** 	Routines for 2D paths (sequences of line segments, also
 **		called `polylines').
 **
 **		This is not a general package for geometric paths, 
 **		which of course include polygons; the emphasis here
 **		is on (for example) usefulness in wire layout.
 **
 ***********************************************************************/

#define	PATHALLOCSIZE(N) \
    (long) ((unsigned) (sizeof(PATH) + \
			(((N)-1) > 0 ? ((N)-1) : 0) \
			* sizeof(Point)))

/*----------------------------------------------------------
 *  String to path / path to string conversion.
 *	External format: 
 *		"(closed, npts, xcoord, ycoord,... )"
 *---------------------------------------------------------*/

PATH *path_in(char *str)
{
    double	coord;
    long	field[2];
    char	*s;
    int	ct, i;
    PATH	*result;
    long	pathsize;
    
    if (str == NULL)
	elog(WARN, "Bad (null) path external representation");
    
    /* read the path header information */
    for (i = 0, s = str; *s && i < 2 && *s != RDELIM; ++s)
	if (*s == DELIM || (*s == LDELIM && !i))
	    field[i++] = atol(s + 1);
    if (i < 1)
	elog(WARN, "Bad path external representation '%s'", str);
    pathsize = PATHALLOCSIZE(field[1]);
    result = (PATH *)palloc(pathsize);
    result->length = pathsize;
    result->closed = field[0];
    result->npts =  field[1];
    
    /* read the path points */
    
    ct = result->npts * 2;	/* two coords for every point */
    for (i = 0;
	 *s && i < ct && *s != RDELIM; 
	 ++s) {
	if (*s == ',') {
	    coord = atof(s + 1);
	    if (i % 2)
		(result->p[i/2]).y = coord;
	    else
		(result->p[i/2]).x = coord;
	    ++i;
	}
    }
    if (i % 2 || i < --ct) {
	PFREE(result);
	elog(WARN, "Bad path external representation '%s'", str);
    } 
    
    return(result);
}


char *path_out(PATH *path)
{
    char 		buf[BUFSIZ + 20000], *result, *s;
    int		i;
    char	tmp[64];
    
    if (path == NULL)
	return(NULL);
    (void) sprintf(buf,"%c%d,%d", LDELIM, 
		   path->closed, path->npts);
    s = buf + strlen(buf);
    for (i = 0; i < path->npts; ++i) {
	(void) sprintf(tmp, ",%G,%G", 
		       path->p[i].x, path->p[i].y);
	(void) strcpy(s, tmp);
	s += strlen(tmp);
    }
    *s++ = RDELIM;
    *s = '\0';
    result = (char *)PALLOC(strlen(buf) + 1);
    (void) strcpy(result, buf);
    
    return(result);
}


/*----------------------------------------------------------
 *  Relational operators.
 *	These are based on the path cardinality, 
 *	as stupid as that sounds.
 *
 *	Better relops and access methods coming soon.
 *---------------------------------------------------------*/

long path_n_lt(PATH *p1, PATH *p2)
{
    return( (p1->npts < p2->npts ) );
}

long path_n_gt(PATH *p1, PATH *p2)
{
    return( (p1->npts > p2->npts ) );
}

long path_n_eq(PATH *p1, PATH *p2)
{
    return( (p1->npts == p2->npts) );
}

long path_n_le(PATH *p1, PATH *p2)
{
    return( (p1->npts <= p2->npts ) );
}

long path_n_ge(PATH *p1, PATH *p2)
{
    return( (p1->npts >= p2->npts ) );
}

/* path_inter -
 *	Does p1 intersect p2 at any point?
 *	Use bounding boxes for a quick (O(n)) check, then do a 
 *	O(n^2) iterative edge check.
 */
long path_inter(PATH *p1, PATH *p2)
{
    BOX	b1, b2;
    int	i, j;
    LSEG seg1, seg2;
    
    b1.xh = b1.yh = b2.xh = b2.yh = (double)DBL_MAX;
    b1.xl = b1.yl = b2.xl = b2.yl = -(double)DBL_MAX;
    for (i = 0; i < p1->npts; ++i) {
	b1.xh = Max(p1->p[i].x, b1.xh);
	b1.yh = Max(p1->p[i].y, b1.yh);
	b1.xl = Min(p1->p[i].x, b1.xl);
	b1.yl = Min(p1->p[i].y, b1.yl);
    }
    for (i = 0; i < p2->npts; ++i) {
	b2.xh = Max(p2->p[i].x, b2.xh);
	b2.yh = Max(p2->p[i].y, b2.yh);
	b2.xl = Min(p2->p[i].x, b2.xl);
	b2.yl = Min(p2->p[i].y, b2.yl);
    }
    if (! box_overlap(&b1, &b2))
	return(0);
    
    /*  pairwise check lseg intersections */
    for (i = 0; i < p1->npts - 1; i++) {
	for (j = 0; j < p2->npts - 1; j++) {
	    statlseg_construct(&seg1, &p1->p[i], &p1->p[i+1]);
	    statlseg_construct(&seg2, &p2->p[j], &p2->p[j+1]);
	    if (lseg_intersect(&seg1, &seg2))
		return(1);
	}
    }
    
    /* if we dropped through, no two segs intersected */
    return(0);
}

/* this essentially does a cartesian product of the lsegs in the
   two paths, and finds the min distance between any two lsegs */
double *path_distance(PATH *p1, PATH *p2)
{
    double *min, *tmp;
    int i,j;
    LSEG seg1, seg2;
    
    statlseg_construct(&seg1, &p1->p[0], &p1->p[1]);
    statlseg_construct(&seg2, &p2->p[0], &p2->p[1]);
    min = lseg_distance(&seg1, &seg2);
    
    for (i = 0; i < p1->npts - 1; i++)
	for (j = 0; j < p2->npts - 1; j++)
	    {
		statlseg_construct(&seg1, &p1->p[i], &p1->p[i+1]);
		statlseg_construct(&seg2, &p2->p[j], &p2->p[j+1]);
		
		if (*min < *(tmp = lseg_distance(&seg1, &seg2)))
		    *min = *tmp;
		PFREE(tmp);
	    }
    
    return(min);
}


/*----------------------------------------------------------
 *  "Arithmetic" operations.
 *---------------------------------------------------------*/

double *path_length(PATH *path)
{
    double	*result;
    int	ct, i;
    
    result = PALLOCTYPE(double);
    ct = path->npts - 1;
    for (i = 0; i < ct; ++i)
	*result += point_dt(&path->p[i], &path->p[i+1]);
    
    return(result);
}



double path_ln(PATH *path)
{
    double	result;
    int	ct, i;
    
    ct = path->npts - 1;
    for (result = i = 0; i < ct; ++i)
	result += point_dt(&path->p[i], &path->p[i+1]);
    
    return(result);
}
/***********************************************************************
 **
 ** 	Routines for 2D points.
 **
 ***********************************************************************/

/*----------------------------------------------------------
 *  String to point, point to string conversion.
 *	External form:	"(x, y)"
 *---------------------------------------------------------*/

Point *point_in(char *str)
{
    char	*coord[POINTNARGS], *p, *r;
    int	i;
    Point	*result;
    
    if (str == NULL)
	elog(WARN, "Bad (null) point external representation");
    
    if ((p = (char *)strchr(str, LDELIM)) == (char *)NULL)
	elog (WARN, "Bad point external representation '%s'",str);
    for (i = 0, p++; *p && i < POINTNARGS-1 && *p != RDELIM; p = r+1)
	if ((r = (char *)strchr(p, DELIM)) == (char *)NULL)
	    elog (WARN, "Bad point external representation '%s'",str);
	else	
	    coord[i++] = p;
    if ((r = (char *)strchr(p, RDELIM)) == (char *)NULL)
	elog (WARN, "Bad point external representation '%s'",str);
    coord[i++] = p;
    
    if (i < POINTNARGS - 1)
	elog(WARN, "Bad point external representation '%s'",str);
    result = PALLOCTYPE(Point);
    result->x = atof(coord[0]);
    result->y = atof(coord[1]);
    return(result);
}

char *point_out(Point *pt)
{
    char	*result;
    
    if (pt == NULL)
	return(NULL);
    result = (char *)PALLOC(40);
    (void) sprintf(result, "(%G,%G)", pt->x, pt->y);
    return(result);
}


Point *point_construct(double x, double y)
{
    Point	*result;
    
    result = PALLOCTYPE(Point);
    result->x = x;
    result->y = y;
    return(result);
}


Point *point_copy(Point *pt)
{
    Point	*result;
    
    result = PALLOCTYPE(Point);
    result->x = pt->x;
    result->y = pt->y;
    return(result);
}


/*----------------------------------------------------------
 *  Relational operators for Points.
 *	Since we do have a sense of coordinates being
 *	"equal" to a given accuracy (point_vert, point_horiz), 
 *	the other ops must preserve that sense.  This means
 *	that results may, strictly speaking, be a lie (unless
 *	EPSILON = 0.0).
 *---------------------------------------------------------*/

long point_left(Point *pt1, Point *pt2)
{
    return( FPlt(pt1->x, pt2->x) );
}

long point_right(Point *pt1, Point *pt2)
{
    return( FPgt(pt1->x, pt2->x) );
}

long point_above(Point *pt1, Point *pt2)
{
    return( FPgt(pt1->y, pt2->y) );
}

long point_below(Point *pt1, Point *pt2)
{
    return( FPlt(pt1->y, pt2->y) );
}

long point_vert(Point *pt1, Point *pt2)
{
    return( FPeq( pt1->x, pt2->x ) );
}

long point_horiz(Point *pt1, Point *pt2)
{
    return( FPeq( pt1->y, pt2->y ) );
}

long point_eq(Point *pt1, Point *pt2)
{
    return( point_horiz(pt1, pt2) && point_vert(pt1, pt2) );
}

/*----------------------------------------------------------
 *  "Arithmetic" operators on points.
 *---------------------------------------------------------*/

long pointdist(Point *p1, Point *p2)
{
    long result;
    
    result = point_dt(p1, p2);
    return(result);
}

double *point_distance(Point *pt1, Point *pt2)
{
    double	*result;
    
    result = PALLOCTYPE(double);
    *result = HYPOT( pt1->x - pt2->x, pt1->y - pt2->y );
    return(result);
}


double point_dt(Point *pt1, Point *pt2)
{
    return( HYPOT( pt1->x - pt2->x, pt1->y - pt2->y ) );
}

double *point_slope(Point *pt1, Point *pt2)
{
    double	*result;
    
    result = PALLOCTYPE(double);
    if (point_vert(pt1, pt2))
	*result = (double)DBL_MAX;
    else
	*result = (pt1->y - pt2->y) / (pt1->x - pt1->x);
    return(result);
}


double point_sl(Point *pt1, Point *pt2)
{
    return(	point_vert(pt1, pt2)
	   ? (double)DBL_MAX
	   : (pt1->y - pt2->y) / (pt1->x - pt2->x) );
}

/***********************************************************************
 **
 ** 	Routines for 2D line segments.
 **
 ***********************************************************************/

/*----------------------------------------------------------
 *  String to lseg, lseg to string conversion.
 *	External form:	"(id, info, x1, y1, x2, y2)"
 *---------------------------------------------------------*/

LSEG *lseg_in(char *str)
{
    char	*coord[LSEGNARGS], *p;
    int	i;
    LSEG	*result;
    
    if (str == NULL)
	elog (WARN," Bad (null) box external representation");
    
    if ((p = (char *)strchr(str, LDELIM)) == (char *)NULL)
	elog (WARN, "Bad lseg external representation '%s'",str);
    for (i = 0, p = str; *p && i < LSEGNARGS && *p != RDELIM; p++)
	if (*p == DELIM || (*p == LDELIM && !i))
	    coord[i++] = p + 1;
    if (i < LSEGNARGS - 1)
	elog (WARN, "Bad lseg external representation '%s'", str);
    result = PALLOCTYPE(LSEG);
    result->p[0].x = atof(coord[0]);
    result->p[0].y = atof(coord[1]);
    result->p[1].x = atof(coord[2]);
    result->p[1].y = atof(coord[3]);
    result->m = point_sl(&result->p[0], &result->p[1]);
    
    return(result);
}


char *lseg_out(LSEG *ls)
{
    char	*result;
    
    if (ls == NULL)
	return(NULL);
    result = (char *)PALLOC(80);
    (void) sprintf(result, "(%G,%G,%G,%G)",
		   ls->p[0].x, ls->p[0].y, ls->p[1].x, ls->p[1].y);
    return(result);
}


/* lseg_construct -
 *	form a LSEG from two Points.
 */
LSEG *lseg_construct(Point *pt1, Point *pt2)
{
    LSEG	*result;
    
    result = PALLOCTYPE(LSEG);
    result->p[0].x = pt1->x;
    result->p[0].y = pt1->y;
    result->p[1].x = pt2->x;
    result->p[1].y = pt2->y;
    result->m = point_sl(pt1, pt2);
    
    return(result);
}

/* like lseg_construct, but assume space already allocated */
void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
{
    lseg->p[0].x = pt1->x;
    lseg->p[0].y = pt1->y;
    lseg->p[1].x = pt2->x;
    lseg->p[1].y = pt2->y;
    lseg->m = point_sl(pt1, pt2);
}

/*----------------------------------------------------------
 *  Relative position routines.
 *---------------------------------------------------------*/

/*
 **  find intersection of the two lines, and see if it falls on 
 **  both segments.
 */
long lseg_intersect(LSEG *l1, LSEG *l2)
{
    LINE *ln;
    Point *interpt;
    long retval;
    
    ln = line_construct_pp(&l2->p[0], &l2->p[1]);
    interpt = interpt_sl(l1, ln);
    
    if (interpt != NULL && on_ps(interpt, l2)) /* interpt on l1 and l2 */
	retval = 1;
    else retval = 0;
    if (interpt != NULL) PFREE(interpt);
    PFREE(ln);
    return(retval);
}

long lseg_parallel(LSEG *l1, LSEG *l2)
{
    return( FPeq(l1->m, l2->m) );
}

long lseg_perp(LSEG *l1, LSEG *l2)
{
    if (! FPzero(l1->m))
	return( FPeq(l2->m / l1->m, -1.0) );
    else if (! FPzero(l2->m))
	return( FPeq(l1->m / l2->m, -1.0) );
    return(0);	/* both 0.0 */
}

long lseg_vertical(LSEG *lseg)
{
    return( FPeq(lseg->p[0].x, lseg->p[1].x) );
}

long lseg_horizontal(LSEG *lseg)
{
    return( FPeq(lseg->p[0].y, lseg->p[1].y) );
}


long lseg_eq(LSEG *l1, LSEG *l2)
{
    return( FPeq(l1->p[0].x, l2->p[0].x) &&
	   FPeq(l1->p[1].y, l2->p[1].y) &&
	   FPeq(l1->p[0].x, l2->p[0].x) &&
	   FPeq(l1->p[1].y, l2->p[1].y) );
}


/*----------------------------------------------------------
 *  Line arithmetic routines.
 *---------------------------------------------------------*/

/* lseg_distance -
 *	If two segments don't intersect, then the closest
 *	point will be from one of the endpoints to the other
 *	segment.
 */
double *lseg_distance(LSEG *l1, LSEG *l2)
{
    double	*d, *result;
    
    result = PALLOCTYPE(double);
    if (lseg_intersect(l1, l2)) {
	*result = 0.0;
	return(result);
    }
    *result = (double)DBL_MAX;
    d = dist_ps(&l1->p[0], l2);
    *result = Min(*result, *d);
    PFREE(d);
    d = dist_ps(&l1->p[1], l2);
    *result = Min(*result, *d);
    PFREE(d);
    d = dist_ps(&l2->p[0], l1);
    *result = Min(*result, *d);
    PFREE(d);
    d = dist_ps(&l2->p[1], l1);
    *result = Min(*result, *d);
    PFREE(d);
    
    return(result);
}

/* distance between l1, l2 */
double lseg_dt(LSEG *l1, LSEG *l2)
{
    double	*d, result;
    
    if (lseg_intersect(l1, l2))
	return(0.0);
    result = (double)DBL_MAX;
    d = dist_ps(&l1->p[0], l2);
    result = Min(result, *d);
    PFREE(d);
    d = dist_ps(&l1->p[1], l2);
    result = Min(result, *d);
    PFREE(d);
    d = dist_ps(&l2->p[0], l1);
    result = Min(result, *d);
    PFREE(d);
    d = dist_ps(&l2->p[1], l1);
    result = Min(result, *d);
    PFREE(d);
    
    return(result);
}


/* lseg_interpt -
 *	Find the intersection point of two segments (if any).
 *	Find the intersection of the appropriate lines; if the 
 *	point is not on a given segment, there is no valid segment
 *	intersection point at all.
 */
Point *lseg_interpt(LSEG *l1, LSEG *l2)
{
    Point	*result;
    LINE	*tmp1, *tmp2;
    
    tmp1 = line_construct_pp(&l1->p[0], &l1->p[1]);
    tmp2 = line_construct_pp(&l2->p[0], &l2->p[1]);
    result = line_interpt(tmp1, tmp2);
    if (result)
	if (! on_ps(result, l1)) {
	    PFREE(result);
	    result = NULL;
	}
    
    PFREE(tmp1);
    PFREE(tmp2);
    return(result);
}

/***********************************************************************
 **
 ** 	Routines for position comparisons of differently-typed
 **		2D objects.
 **
 ***********************************************************************/

#define	ABOVE	1
#define	BELOW	0
#define	UNDEF	-1


/*---------------------------------------------------------------------
 *	dist_
 *		Minimum distance from one object to another.
 *-------------------------------------------------------------------*/

double *dist_pl(Point *pt, LINE *line)
{
    double	*result;
    
    result = PALLOCTYPE(double);
    *result = (line->A * pt->x + line->B * pt->y + line->C) /
	HYPOT(line->A, line->B);
    
    return(result);
}

double *dist_ps(Point *pt, LSEG *lseg)
{
    double m;                       /* slope of perp. */
    LINE *ln;
    double *result, *tmpdist;
    Point *ip;
    
    /* construct a line that's perpendicular.  See if the intersection of
       the two lines is on the line segment. */
    if (lseg->p[1].x == lseg->p[0].x)
	m = 0;
    else if (lseg->p[1].y == lseg->p[0].y) /* slope is infinite */
	m = (double)DBL_MAX;
    else m = (-1) * (lseg->p[1].y - lseg->p[0].y) / 
	(lseg->p[1].x - lseg->p[0].x);
    ln = line_construct_pm(pt, m);
    
    if ((ip = interpt_sl(lseg, ln)) != NULL)
	result = point_distance(pt, ip);
    else  /* intersection is not on line segment, so distance is min
	     of distance from point to an endpoint */
	{
	    result = point_distance(pt, &lseg->p[0]);
	    tmpdist = point_distance(pt, &lseg->p[1]);
	    if (*tmpdist < *result) *result = *tmpdist;
	    PFREE (tmpdist);
	}
    
    if (ip != NULL) PFREE(ip);
    PFREE(ln);
    return (result);
}


/*
 ** Distance from a point to a path 
 */
double *dist_ppth(Point *pt, PATH *path)
{
    double *result;
    double *tmp;
    int i;
    LSEG lseg;
    
    switch (path->npts) {
    case 0:
	result = PALLOCTYPE(double);
	*result = Abs((double) DBL_MAX);	/* +infinity */
	break;
    case 1:
	result = point_distance(pt, &path->p[0]);
	break;
    default:
	/*
	 * the distance from a point to a path is the smallest distance
	 * from the point to any of its constituent segments.
	 */
	Assert(path->npts > 1);
	result = PALLOCTYPE(double);
	for (i = 0; i < path->npts - 1; ++i) {
	    statlseg_construct(&lseg, &path->p[i], &path->p[i+1]);
	    tmp = dist_ps(pt, &lseg);
	    if (i == 0 || *tmp < *result)
		*result = *tmp;
	    PFREE(tmp);
	}
	break;
    }
    return(result);
}

double *dist_pb(Point *pt, BOX *box)
{
    Point	*tmp;
    double	*result;
    
    tmp = close_pb(pt, box);
    result = point_distance(tmp, pt);
    
    PFREE(tmp);
    return(result);
}


double *dist_sl(LSEG *lseg, LINE *line)
{
    double	*result;
    
    if (inter_sl(lseg, line)) {
	result = PALLOCTYPE(double);
	*result = 0.0;
    } else	/* parallel */
	result = dist_pl(&lseg->p[0], line);
    
    return(result);
}


double *dist_sb(LSEG *lseg, BOX *box)
{
    Point	*tmp;
    double	*result;
    
    tmp = close_sb(lseg, box);
    if (tmp == NULL) {
	result = PALLOCTYPE(double);
	*result = 0.0;
    } else {
	result = dist_pb(tmp, box);
	PFREE(tmp);
    }
    
    return(result);
}


double *dist_lb(LINE *line, BOX *box)
{
    Point	*tmp;
    double	*result;
    
    tmp = close_lb(line, box);
    if (tmp == NULL) {
	result = PALLOCTYPE(double);
	*result = 0.0;
    } else {
	result = dist_pb(tmp, box);
	PFREE(tmp);
    }
    
    return(result);
}


/*---------------------------------------------------------------------
 *	interpt_
 *		Intersection point of objects.
 *		We choose to ignore the "point" of intersection between 
 *		  lines and boxes, since there are typically two.
 *-------------------------------------------------------------------*/

Point *interpt_sl(LSEG *lseg, LINE *line)
{
    LINE	*tmp;
    Point	*p;
    
    tmp = line_construct_pp(&lseg->p[0], &lseg->p[1]);
    p = line_interpt(tmp, line);
    if (p)
	if (! on_ps(p, lseg)) {
	    PFREE(p);
	    p = NULL;
	}
    
    PFREE(tmp);
    return(p);
}


/*---------------------------------------------------------------------
 *	close_
 *		Point of closest proximity between objects.
 *-------------------------------------------------------------------*/

/* close_pl - 
 *	The intersection point of a perpendicular of the line 
 *	through the point.
 */
Point *close_pl(Point *pt, LINE *line)
{
    Point	*result;
    LINE	*tmp;
    double	invm;
    
    result = PALLOCTYPE(Point);
    if (FPeq(line->A, -1.0) && FPzero(line->B)) {	/* vertical */
	result->x = line->C;
	result->y = pt->y;
	return(result);
    } else if (FPzero(line->m)) {			/* horizontal */
	result->x = pt->x;
	result->y = line->C;
	return(result);
    }
    /* drop a perpendicular and find the intersection point */
    invm = -1.0 / line->m;
    tmp = line_construct_pm(pt, invm);
    result = line_interpt(tmp, line);
    return(result);
}


/* close_ps - 
 *	Take the closest endpoint if the point is left, right, 
 *	above, or below the segment, otherwise find the intersection
 *	point of the segment and its perpendicular through the point.
 */
Point *close_ps(Point *pt, LSEG *lseg)
{
    Point	*result;
    LINE	*tmp;
    double	invm;
    int	xh, yh;
    
    result = NULL;
    xh = lseg->p[0].x < lseg->p[1].x;
    yh = lseg->p[0].y < lseg->p[1].y;
    if (pt->x < lseg->p[!xh].x)
	result = point_copy(&lseg->p[!xh]);
    else if (pt->x > lseg->p[xh].x)
	result = point_copy(&lseg->p[xh]);
    else if (pt->y < lseg->p[!yh].y)
	result = point_copy(&lseg->p[!yh]);
    else if (pt->y > lseg->p[yh].y)
	result = point_copy(&lseg->p[yh]);
    if (result)
	return(result);
    if (FPeq(lseg->p[0].x, lseg->p[1].x)) {	/* vertical */
	result->x = lseg->p[0].x;
	result->y = pt->y;
	return(result);
    } else if (FPzero(lseg->m)) {			/* horizontal */
	result->x = pt->x;
	result->y = lseg->p[0].y;
	return(result);
    }
    
    invm = -1.0 / lseg->m;
    tmp = line_construct_pm(pt, invm);
    result = interpt_sl(lseg, tmp);
    return(result);
}

Point *close_pb(Point *pt, BOX *box)
{
    /* think about this one for a while */
    
    return(NULL);
}

Point *close_sl(LSEG *lseg, LINE *line)
{
    Point	*result;
    double	*d1, *d2;

    result = interpt_sl(lseg, line);
    if (result)
	return(result);
    d1 = dist_pl(&lseg->p[0], line);
    d2 = dist_pl(&lseg->p[1], line);
    if (d1 < d2)
	result = point_copy(&lseg->p[0]);
    else
	result = point_copy(&lseg->p[1]);
    
    PFREE(d1);
    PFREE(d2);
    return(result);
}

Point *close_sb(LSEG *lseg, BOX *box)
{
    /* think about this one for a while */
    
    return(NULL);
}

Point *close_lb(LINE *line, BOX *box)
{
    /* think about this one for a while */
    
    return(NULL);
}

/*---------------------------------------------------------------------
 *	on_
 *		Whether one object lies completely within another.
 *-------------------------------------------------------------------*/

/* on_pl -
 *	Does the point satisfy the equation? 
 */
long on_pl(Point *pt, LINE *line)
{
    return( FPzero(line->A * pt->x + line->B * pt->y + line->C) );
}


/* on_ps -
 *	Determine colinearity by detecting a triangle inequality.
 */
long on_ps(Point *pt, LSEG *lseg)
{
    return( FPeq (point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
            point_dt(&lseg->p[0], &lseg->p[1])) );
}

long on_pb(Point *pt, BOX *box)
{
    return( pt->x <= box->xh && pt->x >= box->xl &&
	   pt->y <= box->yh && pt->y >= box->yl );
}

/* on_ppath - 
 *	Whether a point lies within (on) a polyline.
 *	If open, we have to (groan) check each segment.
 *	If closed, we use the old O(n) ray method for point-in-polygon.
 *		The ray is horizontal, from pt out to the right.
 *		Each segment that crosses the ray counts as an 
 *		intersection; note that an endpoint or edge may touch 
 *		but not cross.
 *		(we can do p-in-p in lg(n), but it takes preprocessing)
 */
#define NEXT(A)	((A+1) % path->npts)	/* cyclic "i+1" */

long on_ppath(Point *pt, PATH *path)
{
    int	above, next,	/* is the seg above the ray? */
    inter,		/* # of times path crosses ray */
    hi,		/* index inc of higher seg (0,1) */
    i, n;
    double a, b, x, 
    yh, yl, xh, xl;
    
    if (! path->closed) {		/*-- OPEN --*/
	n = path->npts - 1;
	a = point_dt(pt, &path->p[0]);
	for (i = 0; i < n; i++) {
	    b = point_dt(pt, &path->p[i+1]);
	    if (FPeq(a+b,
		     point_dt(&path->p[i], &path->p[i+1])))
		return(1);
	    a = b;
	}
	return(0);
    }
    
    inter = 0;			/*-- CLOSED --*/
    above = FPgt(path->p[0].y, pt->y) ? ABOVE : 
	FPlt(path->p[0].y, pt->y) ? BELOW : UNDEF;
    
    for (i = 0; i < path->npts; i++) {
	hi = path->p[i].y < path->p[NEXT(i)].y;
	/* must take care of wrap around to original vertex for closed paths */
	yh = (i+hi < path->npts) ? path->p[i+hi].y : path->p[0].y;
	yl = (i+!hi < path->npts) ? path->p[i+!hi].y : path->p[0].y;
	hi = path->p[i].x < path->p[NEXT(i)].x;
	xh = (i+hi < path->npts) ? path->p[i+hi].x : path->p[0].x;
	xl = (i+!hi < path->npts) ? path->p[i+!hi].x : path->p[0].x;
	/* skip seg if it doesn't touch the ray */
	
	if (FPeq(yh, yl))	/* horizontal seg? */
	    if (FPge(pt->x, xl) && FPle(pt->x, xh) &&
		FPeq(pt->y, yh))
		return(1);	/* pt lies on seg */
	    else
		continue;	/* skip other hz segs */
	if (FPlt(yh, pt->y) ||	/* pt is strictly below seg */
	    FPgt(yl, pt->y))	/* strictly above */
	    continue;
	
	/* seg touches the ray, find out where */
	
	x = FPeq(xh, xl)	/* vertical seg? */
	    ? path->p[i].x	
		: (pt->y - path->p[i].y) / 
		    point_sl(&path->p[i],
			     &path->p[NEXT(i)]) +
				 path->p[i].x;
	if (FPeq(x, pt->x))	/* pt lies on this seg */
	    return(1);
	
	/* does the seg actually cross the ray? */
	
	next = FPgt(path->p[NEXT(i)].y, pt->y) ? ABOVE : 
	    FPlt(path->p[NEXT(i)].y, pt->y) ? BELOW : above;
	inter += FPge(x, pt->x) && next != above;
	above = next;
    }
    return(	above == UNDEF || 	/* path is horizontal */
	   inter % 2);		/* odd # of intersections */
}


long on_sl(LSEG *lseg, LINE *line)
{
    return( on_pl(&lseg->p[0], line) && on_pl(&lseg->p[1], line) );
}

long on_sb(LSEG *lseg, BOX *box)
{
    return( on_pb(&lseg->p[0], box) && on_pb(&lseg->p[1], box) );
}

/*---------------------------------------------------------------------
 *	inter_
 *		Whether one object intersects another.
 *-------------------------------------------------------------------*/

long inter_sl(LSEG *lseg, LINE *line)
{
    Point	*tmp;

    tmp = interpt_sl(lseg, line);
    if (tmp) {
	PFREE(tmp);
	return(1);
    }
    return(0);
}

long inter_sb(LSEG *lseg, BOX *box)
{
    return(0);
}

long inter_lb(LINE *line, BOX *box)
{
    return(0);
}

/*------------------------------------------------------------------
 * The following routines define a data type and operator class for
 * POLYGONS .... Part of which (the polygon's bounding box is built on 
 * top of the BOX data type.
 *
 * make_bound_box - create the bounding box for the input polygon
 *------------------------------------------------------------------*/

/* Maximum number of output digits printed */
#define P_MAXDIG 12

/*---------------------------------------------------------------------
 * Make the smallest bounding box for the given polygon.
 *---------------------------------------------------------------------*/
void make_bound_box(POLYGON *poly)
{
    double x1,y1,x2,y2;
    int npts = poly->npts;
    
    if (npts > 0) {
	x1 = poly_min((double *)poly->pts, npts);
	x2 = poly_max((double *)poly->pts, npts);
	y1 = poly_min(((double *)poly->pts)+npts, npts),
	y2 = poly_max(((double *)poly->pts)+npts, npts);
	box_fill(&(poly->boundbox), x1, x2, y1, y2); 
    }
}

/*------------------------------------------------------------------
 * polygon_in - read in the polygon from a string specification
 *              the string is of the form "(f8,f8,f8,f8,...,f8)"
 *------------------------------------------------------------------*/
POLYGON *poly_in(char *s)
{
    POLYGON *poly;
    long points;
    double *xp, *yp, strtod();
    int i, size;
    
    if((points = poly_pt_count(s, ',')) < 0)
	elog(WARN, "Bad polygon external representation '%s'", s);
    
    size = offsetof(POLYGON, pts[0]) + 2 * sizeof(double) * points;
    poly = (POLYGON *) PALLOC(size);
    memset((char *) poly, 0, size);	/* zero any holes */
    
    if (!PointerIsValid(poly))
	elog(WARN, "Memory allocation failed, can't input polygon");
    
    poly->npts = points;
    poly->size = size;
    
    /* Store all x coords followed by all y coords */
    xp = (double *) &(poly->pts[0]);
    yp = (double *) (poly->pts + points*sizeof(double));
    
    s++;				/* skip LDELIM */
    
    for (i=0; i<points; i++,xp++,yp++)
	{
	    *xp = strtod(s, &s);
	    s++;					/* skip delimiter */
	    *yp = strtod(s, &s);
	    s++;					/* skip delimiter */
	}
    make_bound_box(poly);
    return (poly);
}

/*-------------------------------------------------------------
 * poly_pt_count - count the number of points specified in the
 *                 polygon.
 *-------------------------------------------------------------*/
long poly_pt_count(char *s, char delim)
{
    long total = 0;
    
    if (*s++ != LDELIM)		/* no left delimeter */
	return (long) -1;
    
    while (*s && (*s != RDELIM))
	{
	    while (*s && (*s != delim))
		s++;
	    total++;	/* found one */
	    if (*s)
		s++;	/* bump s past the delimiter */
	}
    
    /* if there was no right delimiter OR an odd number of points */
    
    if ((*(s-1) != RDELIM) || ((total%2) != 0))
	return (long) -1;
    
    return (total/2);
}

/*---------------------------------------------------------------
 * poly_out - convert internal POLYGON representation to the 
 *            character string format "(f8,f8,f8,f8,...f8)"
 *---------------------------------------------------------------*/
char *poly_out(POLYGON *poly)
{
    int i;
    double *xp, *yp;
    char *output, *outptr;
    
    /*-----------------------------------------------------
     * Get enough space for "(f8,f8,f8,f8,...,f8)"
     * which P_MAXDIG+1 for each coordinate plus 2
     * for parens and 1 for the null
     *-----------------------------------------------------*/
    output = (char *)PALLOC(2*(P_MAXDIG+1)*poly->npts + 3);
    outptr = output;
    
    if (!output)
	elog(WARN, "Memory allocation failed, can't output polygon");
    
    *outptr++ = LDELIM;
    
    xp = (double *) poly->pts;
    yp = (double *) (poly->pts + (poly->npts * sizeof(double)));
    
    sprintf(outptr, "%*g,%*g", P_MAXDIG, *xp++, P_MAXDIG, *yp++);
    outptr += (2*P_MAXDIG + 1);
    
    for (i=1; i<poly->npts; i++,xp++,yp++)
	{
	    sprintf(outptr, ",%*g,%*g", P_MAXDIG, *xp, P_MAXDIG, *yp);
	    outptr += 2*(P_MAXDIG + 1);
	}
    *outptr++ = RDELIM;
    *outptr = '\0';
    return (output);
}

/*-------------------------------------------------------
 * Find the largest coordinate out of n coordinates
 *-------------------------------------------------------*/
double poly_max(double *coords, int ncoords)
{
    double max;
    
    max = *coords++;
    ncoords--;
    while (ncoords--)
	{
	    if (*coords > max)
		max = *coords;
	    coords++;
	}
    return max;
}

/*-------------------------------------------------------
 * Find the smallest coordinate out of n coordinates
 *-------------------------------------------------------*/
double poly_min(double *coords, int ncoords)
{
    double min;
    
    min = *coords++;
    ncoords--;
    while (ncoords--)
	{
	    if (*coords < min)
		min = *coords;
	    coords++;
	}
    return min;
}

/*-------------------------------------------------------
 * Is polygon A strictly left of polygon B? i.e. is
 * the right most point of A left of the left most point
 * of B?
 *-------------------------------------------------------*/
long poly_left(POLYGON *polya, POLYGON *polyb)
{
    double right, left;
    
    if (polya->npts > 0)
	right = poly_max((double *)polya->pts, polya->npts);
    else
	right = polya->boundbox.xh;
    if (polyb->npts > 0)
	left = poly_min((double *)polyb->pts, polyb->npts);
    else
	left = polyb->boundbox.xl;
    
    return (right < left);
}

/*-------------------------------------------------------
 * Is polygon A overlapping or left of polygon B? i.e. is
 * the left most point of A left of the right most point
 * of B?
 *-------------------------------------------------------*/
long poly_overleft(POLYGON *polya, POLYGON *polyb)
{
    double left, right;
    
    if (polya->npts > 0)
	left = poly_min((double *)polya->pts, polya->npts);
    else
	left = polya->boundbox.xl;
    if (polyb->npts > 0)
	right = poly_max((double *)polyb->pts, polyb->npts);
    else
	right = polyb->boundbox.xh;
    
    return (left <= right);
}

/*-------------------------------------------------------
 * Is polygon A strictly right of polygon B? i.e. is
 * the left most point of A right of the right most point
 * of B?
 *-------------------------------------------------------*/
long poly_right(POLYGON *polya, POLYGON *polyb)
{
    double right, left;
    
    if (polya->npts > 0)
	left = poly_min((double *)polya->pts, polya->npts);
    else
	left = polya->boundbox.xl;
    if (polyb->npts > 0)
	right = poly_max((double *)polyb->pts, polyb->npts);
    else
	right = polyb->boundbox.xh;
    
    return (left > right);
}

/*-------------------------------------------------------
 * Is polygon A overlapping or right of polygon B? i.e. is
 * the right most point of A right of the left most point
 * of B?
 *-------------------------------------------------------*/
long poly_overright(POLYGON *polya, POLYGON *polyb)
{
    double right, left;
    
    if (polya->npts > 0)
	right = poly_max((double *)polya->pts, polya->npts);
    else
	right = polya->boundbox.xh;
    if (polyb->npts > 0)
	left = poly_min((double *)polyb->pts, polyb->npts);
    else
	left = polyb->boundbox.xl;
    
    return (right > left);
}

/*-------------------------------------------------------
 * Is polygon A the same as polygon B? i.e. are all the
 * points the same?
 *-------------------------------------------------------*/
long poly_same(POLYGON *polya, POLYGON *polyb)
{
    int i;
    double *axp, *bxp; /* point to x coordinates for a and b */
    
    if (polya->npts != polyb->npts)
	return 0;
    
    axp = (double *)polya->pts;
    bxp = (double *)polyb->pts;
    
    for (i=0; i<polya->npts; axp++, bxp++, i++)
	{
	    if (*axp != *bxp)
		return 0;
	}
    return 1;
}

/*-----------------------------------------------------------------
 * Determine if polygon A overlaps polygon B by determining if
 * their bounding boxes overlap.
 *-----------------------------------------------------------------*/
long poly_overlap(POLYGON *polya, POLYGON *polyb)
{
    return box_overlap(&(polya->boundbox), &(polyb->boundbox));
}

/*-----------------------------------------------------------------
 * Determine if polygon A contains polygon B by determining if A's
 * bounding box contains B's bounding box.
 *-----------------------------------------------------------------*/
long poly_contain(POLYGON *polya, POLYGON *polyb)
{
    return box_contain(&(polya->boundbox), &(polyb->boundbox));
}

/*-----------------------------------------------------------------
 * Determine if polygon A is contained by polygon B by determining 
 * if A's bounding box is contained by B's bounding box.
 *-----------------------------------------------------------------*/
long poly_contained(POLYGON *polya, POLYGON *polyb)
{
    return box_contained(&(polya->boundbox), &(polyb->boundbox));
}
